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Theories of modified gravity, in both the linear and fully non-linear regime, are often studied under 
the assumption that the evolution of the new (often scalar) degree of freedom present in the theory 
is quasi-static. This approximation significantly simplifies the study of the theory, and one often 
has good reason to believe that it should hold. Nevertheless it is a crucial assumption that should 
be explicitly checked whenever possible. In this paper we do so for the Vainshtein mechanism. 
By solving for the full spatial and time evolution of the Dvali-Gabadadze-Porrati and the Cubic 
Galileon model, in a spherical symmetric spacetime, we are able to demonstrate that the Vainshtein 
solution is a stable attractor and forms no matter what initial conditions we take for the scalar 
field. Furthermore,the quasi-static approximation is also found to be a very good approximation 
whenever it exists. For the best-fit Cubic Galileon model, however, we find that for deep voids at 
late times, the numerical solution blows up at the same time as the quasi-static solution ceases to 
exist. We argue that this phenomenon is a true instability of the model. 


I. INTRODUCTION 

There is overwhelming evidence that the Universe is 
undergoing accelerated expansion. A possible, popular, 
explanation is the existence of a cosmological constant A. 
The A Cold Dark Matter (ACDM) model is, thus far in 
perfect agreement with observations. An alternative is 
that new dynamical degrees of freedom must be at play 
in our Universe, in the form of a new fluid that has been 
dubbed dark energy [T[. Yet another alternative is that 
gravity is modified on cosmological scales [5] ■ 

General Relativity has been exquisitely tested in the 
solar system [3]. Any theory that attempts to modify 
gravity must conform to the stringent limits coming from 
these experiments; this usually requires the presence of 
a screening mechanism mu, i.e. a way of dynamically 
hiding modifications of gravity in regions where we have 
tested general relativity to great precision. The Vain¬ 
shtein mechanism [8] is one such mechanism and relies on 
non-linear differential operators to screen the fifth-force 
in high density regions of spacetime. This mechanism is 
responsible for, the viability of popular modified gravity 
models such as massive gravity [3], DGP na and the 
Galileon [TTHT3] . 

The study of modified gravity models often require 
us to solve complicated, non-linear, partial differential 
equations for the new degrees of freedom. For scales 
that are small compared to the horizon one can often 
apply what is known as the quasi-static approximation 
to significantly simplify the analysis. This approxima¬ 
tion roughly means that we neglect most terms involving 
time-derivatives of perturbed quantities relative to those 
containing spatial derivatives. On linear, sub-horizon, 
scales it has been shown that, for a wide class of mod¬ 
els, we can safely assume the quasi-static approximation. 
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However as was shown in m its is possible to have large 
deviations on scales close to the horizon. In general we 
have that the quasi-static approximation breaks down for 
scales larger than the sound horizon of the dynamical de- 
gree(s) of freedom in question [ITS] . 

In the non-linear regime of structure formation the 
main tool to obtain accurate predictions is A-body sim¬ 
ulations. Over the last decade there have been a number 
of studies of screening mechanisms [T51U5j in the non¬ 
linear regime by performing such simulations. Most of 
these simulations have assumed the quasi-static approx¬ 
imation. Recently a handful of papers have appeared 
where scalar fields have been allowed to evolve dynami¬ 
cally |26fi30] and treated beyond the quasi-static approx¬ 
imation. These studies have demonstrated that, in most 
situations of cosmological interest, the approximation is 
very good. Nevertheless interesting non-static effect have 
also been found in certain models, such as, for example, 
the formation of domain walls m- 

For the Vainshtein mechanism, which is the topic 
of this paper, the quasi-static approximation has been 
shown to work well for linear perturbations (323, but for 
the non-linear regime the only study we are aware of 
looking at non-static effects was in |22] for the Dvali- 
Gabadadze-Porratti (DGP) model. It was shown that the 
quasi-static approximation is a self-consistent approxi¬ 
mation in the sense that by assuming it and solving the 
approximate equations, one could then calculate the ef¬ 
fect of the neglected terms and show that they were in¬ 
deed negligible. 

Another interesting class of models with Vainshtein 
screening are the so-called Galileon models (see e.g. (311 - 
00]). The Cubic Galileon model, being the simplest 
Galileon, has been found to be a viable alternative to 
ACDM for explaining cosmic accelerator^ (411 1 142 ]. How¬ 
ever, the best-fit Cubic Galileon model has a peculiar 


1 There are some issues related to the ISW effect in the Galileon 
models that might be is in tension with observations. 
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property: if we study the theory in the quasi-static ap¬ 
proximation for a spherical symmetric body in a Fried¬ 
man Robertson Walker metric then we find that, for some 
configurations (depending on the model parameters and 
the object in question), the analytic field-profile ceases 
to exist. The same effect has been found in iV-body sim¬ 
ulations of this model. It has been speculated m that 
this might be an artifact of the quasi-static approxima¬ 
tion and that if non-static effects are taken into account 
then this problem could be solved. This is one of the 
questions we aim to answer here. 

In this paper we present spherical symmetric cosmo¬ 
logical simulations of the DGP and the Cubic Galileon 
model beyond the quasi-static approximation. We want 
to answer questions like: does the Vainshtein mechanism 
always form? Is the quasi-static approximation always 
valid in the non-linear regime? Are other approximations 
widely used in simulations of such models, like neglect¬ 
ing terms with metric potentials whenever they do not 
have a V 2 in front of them, also valid? In all the ques¬ 
tions above we naively expect the answer to be yes. Our 
aim is to put these questions to rest and place the quasi¬ 
static approximation on a firm footing by demonstrating 
and quantifying its validity in the non-linear regime. 

The setup of this paper is as follows. In Sec. © 
we intr oduce the DGP and the Cubic Galileon model, 
in Sec. (IIC) we discuss the Vainshtein mechanism and 
the quasi-static approximation. In Sec. (IllI we present 
the equations needed to go beyond the quasi-static ap¬ 
proximation and in Sec. (IV ) we discuss the numerical 
implementation before presenting code-tests in Sec. fv| ) 
and finally the results of our simulations in presented in 
Sec. 

Unless stated otherwise we always work in units of 
c = h = 1 and Mpi = \Jhc/ ( 87 rG) = 1. Unless stated 
otherwise we use the metric sign-convention (—,+,+,+). 


II. THE VAINSHTEIN MECHANISM 

The Vainshtein Mechanism is the name of the mecha¬ 
nism with which the Galileon and DGP model is able to 
evade the stringent local gravity constraints in the solar- 
system. For a comprehensive review of how the mecha¬ 
nism (and other screening mechanisms) works in detail 
see [T . We will focus on how it arises in the cosmological 
context and in particular its relation with the quasi-static 
limit. In this section we review two simple models that 
have the Vainshtein mechanism. 


A. The DGP model 

The DGP model [10] is a modified gravity model where 
we are confined to live in a four-dimensional brane, em¬ 
bedded in a five-dimensional spacetime. The gravita¬ 


tional action is given by 

f , f r R( 5 ) 

s DG p=y<ix^—+y dI ^_ w , (!) 

where and denote quantities on the brane and in 
the bulk respectively. The relative sizes of the two gravi¬ 
tational constants is a parameter of the model known as 
the crossover scale, r c , 


1G(S) 
r ° ~ 2 ~G~ 


( 2 ) 


For scales smaller than r c gravity looks four-dimensional. 
This model, in its simplest form, has two branches of solu¬ 
tions. The first branch, the normal branch, requires dark 
energy to produce cosmic acceleration of the Universe. 
The second branch, the self-accelerating branch, is able 
to produce acceleration without a cosmological constant 
but is in tension with observations |43j and is plagued by 
theoretical problems associated with the propagation of 
ghosts [35]. For our purposes, i.e. to test the validity of 
the quasi-static approximation, these problems are not 
important and we will look at both branches. 

Assuming a homogeneous and isotropic matter distri¬ 
bution, an empty Minkowski bulk and a (spatially) flat 
brane leads to the Friedmann equation 


A 2 ± ^ = 
rv 


87 tG 

~T~ Pm 


(3) 


where the sign ± refers to the two branches mentioned 
above. In the accelerating branch we obtain 


^ — \/+ V^rc + 3 (4) 

riQ 

where f2 rc = ir 2 H -i and = 1 — r Ho . Looking at 
perturbations about the Friedmann solution, an addi¬ 
tional degree of freedom appears which is associated with 
displacement of the brane. This is the so-called brane- 
bending mode 0 and is determined by the field equation 

|33H3B] 


r 2 RttC 

□</> + [(n0) 2 - (v M v,y») 2 ] = —s P (5) 

where 


/?(a) 


1 ± 217(a) 



j 

3 H 2 (a) J 


( 6 ) 


Since 0 is a perturbative quantity the background value 
of this is simply 0 = 0. The gravitational potential T in 
the model is the sum of the standard Newtonian potential 
'Fjv and the brane-bending mode, 


4 / = 4/at + ^ (7) 

giving rise to a fifth-force on matter-particles F 0 = | V 0 . 

In the case of the normal branch (+ sign in /3) we will in 
this paper simply assume that the background evolution 
of the model is the same as in ACDM. 
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B. The Cublic Galileon 


The Cublic Galileon model also containis a Vainshtein 
mechanism and the resulting equations are very similar 
to that of DGP. However there are some important dif¬ 
ferences as we shall see below. The action is given bj 0 


S= d^Xyf^g 


R Co _ C 3 „ 

—2" £2 — 2 ^ 3 ~ ^ matter 


where 02,03 are dimensionless parameters, 
£2 = (V <^) 2 


£3 = M3 (V ^ 


( 8 ) 

(9) 

( 10 ) 


and M 3 = My>\Hq. The Einstein equations following 
from the action are 


1 _ rr imatter 1 rp(p{2) 1 /-t-i<^>( 3 ) 

'^ r l±v — 1 LLV ' 1 fiu ' 1 nv 


( 11 ) 


where the energy-momentum tensor of the scalar field is 
given by 


=c 2 
2 c : 


^3^(V</>) 2 


( 12 ) 




M 3 

— 2 Va0V( ai </>V„)V > 


(13) 


A variation of the action with respect to <f> gives us a 
modified Klein-Gordon equation: 

+ ^3 [( D </>) 2 - ( V , V , 0) 2 - i ^ V ^] = 0 

(14) 

In a flat Friedmann-Robertson-Lemaitre-Walker metric 
(FRLW) the Einstein equations gives us the Friedmann 
equations 


3 H 2 = p m + + 6 c 3 


1 -2 


H(j) 


H + H = —- 


Pm + 2c 2 (t) + 6 c 3 —3 (#</>-</>) 

n o 


(15) 

(16) 


and the equation of motion reduces to 

6 C „ *2 ^_2 . 

4 i + 377</> + —^ (3£r 2 <(. + 2(j>(j)H + <j> H) = 0 (17) 

c 2 U 0 


These equations have a late time de Sitter attractor given 
by H(f> = To obtain the correct dark energy 

density today we need 


— 


C2</>0 

6 HI 


2 c 3 


U 3 
^0 


(18) 


and by evaluating this expression using the attractor we 

c 3 

find = 1 — = — pp. The value of c 2 can be fixed 

by perfoming a rescaling of </> so we can without loss of 
generality choos^]c 2 = — 1. We then choose 


C3 = 


6^6(1-fi m ) 


(19) 


so the model has the same number of free parameters as 
ACDM. In the following we will use the best-fit param¬ 
eters from the analysis in m which can be translated 
into c 2 = —1 and c 3 = 0.08 corresponding to = 0.277 
and = 0.723. 

In Sec. (IllI we will discuss the specific form of the 


equations above when we go to a spherical symmetric 
spacetime. 


C. The Vainshtein Mechanism in the Quasi-Static 
Approximation 

The quasi-static approximation can be defined as the 
limit of the theory for which time derivatives of the per¬ 
turbations around a given background can be neglected. 

Roughly this reduces to making the replacement </> —» 4> 
where an overline denotes a background quantity. For 
cosmological perturbation theory we also neglect terms 
proportional to $ where $ is any metric potential. The 
discussion below will be for the Cubic Galileon, but the 
same equations with some small modifications, also apply 
for the DGP model. 

Applying the quasi-static approximation to the Klein 
Gordon equation Eq. (14) for the Cubic Galileon we find 


V 2 d> + 


1 


3/3ii7 0 2 a 2 


[(V 2 0 ) 2 - (V 4 V^) 2 ] = —Spr, 


where 


Pi( a ) = 


6 c 3 


j_4 


— C 2 — Acs^cj) “K 2,H(j)') 2 Cr$(f> 


&(a) = 


( 20 ) 

( 21 ) 

( 22 ) 


2 Note that in this subsection, and only in this subsection, we use 
the metric sign-convention (+, to be consistent with the 

convention used in m- 


3 The sign here is not determined by the rescaling. We choose ’ 
which is the same as in the best-fit cosmological model found in 

m- 
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are functions of the scale factor. The DGP model fits also 
into this description by taking /3i(a) = /3(a)/ {t c Hq) 2 and 

02 ( a ) = 0 (a)- 

The above equations are general, but if we assume 
spherical symmetry then we can simplify further as it 
is integrable; integrating it over r we obtain 

2 a 2 5M(r) 

12tt/ 3 2 r 3 1 J 

where 5M(r) = 47T 5p m r 2 dr. This is an algebraic 

equation in <p'/r which can be solved to give us 


r ZPiHZa 2 \ r 



- 1 


(24) 


where ry = is the Vainshtein radius of an object 

with Schwarchild radiutj^Jrg = 2 GSM. This equation can 
be rewritten as 


2r 3 


= a x 


with a = 


1 + ^-1 


GSM(r ) 


(25) 


If we now take, see Eq. (34), V 2 4/ ~ V 2 \I/jv — ^fls-V 2 ^, 
1 ^0 

we have that the gravitational force F = V’F in the 
Newtonian limit predicted by the Galileon model can be 
rewritten as F = Fn + F^ = ygFFjv where Fn is the 
Newtonian expression for the force and where the effec¬ 
tive gravitational constant, which is a function of both 
time and scale, satisfies 


G e ff 


= 1 - 


C 3 </> 2r 3 

3^7 x 


i 


r 3 

J- 1 


(26) 


The equivalent expression for the DGP model can be 

_!_2 

recovered by making the replacement —> 3 ^ a ^ . 

Outside the Vainshtein radius r ry gravity is modified 
as 

Geff-G 


G 


~ <7 


(27) 


where a = for the Cubic Galileon and er = 

for the DGP model. Within the Vainshtein radius on the 
other hand (r <C ry) we have 

Geff-G 


G 


~ cr X 


r \ 3/2 

- < CT 

r v J 


(28) 


If an object it massive enough as to have a large Vain¬ 
shtein radius then, for r < ry General Relativity is re¬ 
covered. 


D. Breakdown of the quasi-static Vainshtein 
solution in the Cubic Galileon model 


The quasi-static spherical symmetric solution derived 
in the previous subsection is only valid when 


1 + 


'v 


3 >° 


(29) 


otherwise we get a complex solution for <//. This condi¬ 
tion translates into 


9010 2<t 3 3 fg S„ 




> — 


r 2 dr 


= - (Sm) 


(30) 


For overdensities this condition is always satisfice^ as 

S m > 0 and /3 i/3 2 oc (^/3i/<frj > 0. For DGP we also 

have 0i02 = 0 2 (a)/(r c H o ) 2 > 0. Furthermore, in both 
the self-accelerating and the normal branch DGP model 
we have that the critical value for the density contrast 


8 C < — | < — 1 at all times, see Fig. ([2|, so the quasi¬ 
static solution always exists. 

For the Cubic Galileon model, however, the condition 
breaks down at late times in voids. For a void with den¬ 
sity contrast <5 q the condition reads 


^0 > — 


9/3i/3 2 a 3 

8fl m 


(31) 


In Fig. ([lj we show the time evolution of this condition 
as a function of the scale-factor. At a > 0.8 the condition 
starts to break down for the deepest voids <5 0 — — 1 and 
at the present time the condition is violated for all voids 
with S 0 < —0.5. 

This was first noted in pi] where the IV-body simula¬ 


tions, based on solving Eq. (201, broke down in deep cos¬ 


mic voids close to the present time. The authors choose 
to make an ad hoc fix to be able to integrate the equa¬ 
tions until today with the argument that even though it 
breaks down in voids, whatever happens there should not 
significantly alter clustering statistics such as the halo 
mass function and the matter power-spectrum. It was 
also pointed out that this effect could to be due to the 
assumptions made, namely the quasi-static approxima¬ 
tion. 

We would like to know what happens to the solution 
when we enter this regime. Does it blow up or do the 
terms which were neglected when going to the quasi¬ 
static limit kick in and save the day? To answer such 
question we must move beyond the quasi-static limit 
and look at solutions to the equations when the time- 
evolution is taken properly into account. This will be 
the subject in the next section. 


4 The mass SM is measured relative to the cosmic mean so a void¬ 
like structure has a negative rg. 


5 This is only for the Cubic and Quartic Galileon. For the general 
Galileon model this problem can also occur for 0(1) overdensities 
as shown in m ■ 
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FIG. 1. The allowed density contrast S c of a void versus scale 
factor a for which the condition Eq. (311 is satisfied for the 
Cubic Galileon model. For a > 0.81 we start to see 5 C > —1 
and in the far future (a —> oo) we have 5 C —► 0. At the 
present time, a = 1, the analytic quasi-static field profile does 
not exist today for voids with density contrast S < —0.52. 



FIG. 2. The allowed density contrast 5 C o f a void versus scale 
factor a for which the condition Eq. (311 is satisfied for the 


DGP model. 


III. BEYOND THE QUASI-STATIC 
APPROXIMATION 

In this section we will derive the full evolution equa¬ 
tions, including metric perturbations, for our models. We 
are, for simplicity going to focus on spherically symmetric 
systems relevant for cosmology. The presentation in this 
section will focus on the Cubic Galileon. The non-linear 
terms in the Galileon models vanish for ID configurations 
so the spherical symmetric case is the simplest setup for 
which we have the Vainshtein mechanism present. The 
equations and analysis for the DGP model is very similar 
so we will present it in Appendix [B] instead. 

We work with the perturbed FRLW metric in the New¬ 
tonian Gauge 

ds 2 = -dt 2 { 1 + 2\P) + a 2 (f)( 1 - 2$) (dr 2 + r 2 dQ 2 ) (32) 

To begin with we are going to make as few approxima¬ 
tions as possible: 

• The metric potentials <E>, T are expanded to first 
order in perturbation theory. 

• We neglect all factors of $, T that don’t contain 
a derivative, i.e. (1 + <f>) ~ 1. These terms only 
renormalize the terms we already have by a very 
small factor < 1 + 0(1O -5 ), for the cosmological 
case, so it is safe to ignore them. 

• The scalar field (f> = (f> + S(f> is treated exactly, no 
perturbative expansion and no quasi-static approx¬ 
imation is applied. 

In Appendix[C]we present expressions for the ChristofFel- 
symbols, the Riemann tensor and derivative operators on 
the scalar field needed to derive the equations presented 
below. 


A. The dimensionless field equations 


Before we present the equations we will solve we will 
first introduce our dimensionless code-units. Our code 
variables are defined as 


x = log a, 




Q = 


1 d(f> 
H 0 dt 


(33) 


where R is the size of the simulation box at the present 
time and a = , J . Further we will also use V 2 d> = 

ayriQti) y T 

j 2 (y 2 w) = ^ to simplify the notation. An 

overbar will denote a background quantity and E = 
is the Hubble factor in units of the present value. We have 
also implemented physical coordinates y p h ys = ^ = ay 
and to get the equations in these coordinates the only 
change one has to make is the replacement a —> aer in all 
the equations below. 

The first equation that is needed is the perturbed part 
of the Einstein equation, dG^ 00 ^ = 5T ^ 00 - ) . This gives us 
the Poisson equation for T which in code units can be 
written 
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C 2 , 


-^-2 


^3, 


a 2 V 2 * =^ m a- a «5 m + ^(Q 2 -Q‘ + <J 2 co 2 ) + 3 Ec 3 (Q A - Q“) - c 3 (Q 2 - aV)a 2 V 2 ' 


— c 3 EQa 2 uj 2 — 


2 cr 4 u; 3 c 3 


+ [3£ 2 - 3£Q 3 c 3 + c 3 EQo 2 u 2 \ + c 3W a 2 [Q 2 + aV] (34) 


For the anisotropic stress we take 5G r r — l< 5G\ = 5T$ — YT\ which gives us: 


cr 2 V 2 (4> -$)-—(*-$)' = 2c 3 <t 2 w 2 3AQ + (7 2 V 2 </> - 


4cr 2 a; C 2 , „ (dQ 


V 


, E . , 
2c 3 V dx 


— 4c 3 ct 2 CliQ 


+2c 3 ct“w 


— EQuj 


d$ 

dx 


(*S\- 


+ (2Q - a 2 io 2 ) — - 3.EQu; — + 3cr 2 w 


\dy J 


dT\ 


dx J 


dQ 

dy 

d$\ 

dy). 


(35) 


r 


The right hand side vanishes if oj = 0 and, in that case, we have the Klein Gordon equation for (j) which can be 

1 d4> , 

Hr) dt ■ 


we expect <F ~ 'F to hold to a high accuracy. Finally written as an evolution equation for Q = -4-^: 


E 


C 2 


12.EQ + — - 4cr V 2 0 - 12.EQ ( — I + 4cr w ( — ) 


d\F\ 


3EQc2 


c 3 


c 3 


8 Ea 2 u [ - 4 (t 2 ) + cr 2 V 2 (j 


d\F\ 


dQ 


dy 


dy 


dy J \dx 


‘*+8EQ- s ^ 
c 3 y 




- 18£ 2 Q 2 + 2 


dE 

dx 


Ea z cu z + 2 E z a z co z + 


12cr 4 a; 2 

y 2 


- 2 ct 2 V 2 (’F - $) + 2 ( t 2 V 2 'F(Q 2 + crV) 


■£Ui 


d$ 

dx 


,'dtf\ , /d$ 

E - 42 | *) +<, - 4 n* 


where ^4 i_ 6 are functions given in Appendix [A] 


. ,'dtf\ r , ( d 1 V F \ _, 2 . /d 2 \F 

^ ( dyj (dydx) A (dx^ 


r 


(36) 


Removing the perturbations u, $ -> 0, Q -> Q we get the homogenous scalar field equation for the background 


dQ 

dx 


3QEc 2 . a -^2 


-^-2 „ dE 


c 3 


- 18Q E z -6Q E 


dx 


x | | j + 12£Q^ E 


(37) 


which can be written using Eq. (151 as 
3QEc 2 


dQ 

dx 


- 12 Q~E 2 + 3 n m e- 3x Q 2 + 2c 2 Q? + 6c 3 Q 5 E \ x ( ( ^ + 12FiQ + 6c 3 Q 4 1 E 


c 3 


(38) 


This equation can be seen to be identical to Eq. and 
serves as a simple consistency check. Lastly we need the 


Hubble equation Eq. (151 on a dimensionless form: 


7^3 


E = c 3 Q u + \/c 2 Q 6 + -jvQ 2 + f2 m e 3x (39) 
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The equations above form a closed system, as long as 
fim(y,x) is given, and are all we need for our numerical 
implementation. 


B. Conserved Charge 

The Galileon action is invariant under the symmetry 
<j> 4> + c + d^x 11 and consequently, through Noether’s 

theorem, there is a conserved current in terms of 
which which the equation of motion can be written 

V ^ = l ( ^ J ° ) + ^ (v/ ^ jr) = ° (40) 

where the first equality only holds for spherical symme¬ 
try. Associated with this current is a conserved change 
densitj0 


_ f V=gJ°dr 

PNoether — r o , 

J r z ar 


(41) 


where 


J “ = (f + w a4 ) - (42) 


and in particular 


^J°=e x ^y Q (yy + 6EQ - 2ct 2 V 2 c/> 
+2jA<rV - 3 Q 2 ) - 2 EaW + 2 a 2 w^ 

dx ay 

+2Qa 2 u>y~(^ - T) 

dy 


(43) 


For the background the conserved charge density be¬ 
comes 


PNoether ^ 


3 x ( ^2 
C3 


6 EQ Q 


(44) 


C. What terms are expected to be small? 


One of the aims with this paper is to rigorously classify 
what are terms in the (long) equations above that can be 
safely neglected. We will first go through the standard 
arguments for why some terms should be small, then we 
will study this numerically in the next section. 

All the terms in Ai. 2,6 can be found elsewhere in the 
equation of motion so these terms are expected to be 

C 1. In a cosmo¬ 


1 1 

1 d * 1 


1 dx \ 

1 dx \ ’ 

dx 2 


logical context the potentials Y = 4>, T usually evolve 
over a Hubble time which makes ^ ~ y -C 1. The 
terms containing derivatives of the gravitational poten¬ 
tial need a closer inspection. If we define ( 5 m (r )) to be 
the average matter density within a radius r we roughly 
have 


dr 



Ilm (S m (r)) 


r 


(46) 


and likewise for $. This leads to the constraint 


a 6000 
( 8m(r )) 


Mpc /h 


(47) 


which is satisfied for most cosmological applications. We 
also expect perturbations in the scalar field to be small. 
In particular we expect 


\6cj)\ < 1 , 



d(f> 

dt 


5cf) 




(48) 


If the first condition is violated then we are outside the 
realm of the effective theory which we are working with. 
If the second or third condition is violated it means 
that the clustering of the scalar field will back-react and 
thereby changing the background evolution of the Uni¬ 
verse. 

Lastly we have the anisotropic stress. If there is no 
perturbations in the scalar field this vanishes identically 
and in general we have 

V2($-4-)^0(w 2 Q) (49) 


By taking the derivative of this equation, dpN ^° ther = 0, 
we recover the background equation of motion Eq. 

We can use the conservation of the Noether charge den¬ 
sity as a test of our code by monitoring the constancy 
of 

£ _ PNoether ~ PNoether 
PNoether 


Given that w 2 ~ T 2 and EQ = 0(1) we also expect the 
anisotropic stress to be negligible. 

We finish by pointing out that there is one term ne¬ 
glected in the quasi-static limit that could potentially be 

t 2 i 

large and that is the cross term We would gener¬ 

ally expect jbA ~ which is comparable with other 
terms in the equation of motion. 


The integrals here are to be interpreted as first taken over a 
fixed volume in co-moving coordinates and then taking the limit 
of the volume going to infinity. For numerical simulations with 
periodic boundary conditions the integrals are simply taken over 
the simulation volume. 


IV. NUMERICAL IMPLEMENTATION 

The fundamental variables for the scalar field in our 
code are the ’position’ (f> and ’velocity’ Q = j^cj) which 
are discretised on a linearly spaced grid going from y = 0 
to y = 1 with N grid cells. 




















A. Leap-Frog Integrator 

Our first method to solve the system is to propagate 
the variables {(/>, Q} in time using a staggered Leap Frog 
algorithm. Starting with Q n -\/2 and <p n we first propa¬ 
gate Q one step using 


for second order accuracy. We have also implemented 
and tested higher order discretizations, but our second 
order implementation was found to be well behaved and 
no visible difference in the results was found when we 
included higher order terms. 


f)„ a " (50) 

and then use the result to propagate <f> using 

4>n+ 1 = 4>n + [QE\ n+ i/ 2 Ax (51) 

From </>„+! we calculate u n +i = d< ^ +1 and other spa¬ 
tial derivatives depending on fusing a five-point stencil 
which is fourth order accuratq^] This is then again used 
in the next time-step to evaluate . 


Qn+ 1/2 — Qn—1/2 + ^ 


B. Newton Gauss Seidel Integrator 


The Klein Gordon equation Eq. (361 and Eq. (B11 can 
be written on the schematic form £ = A(cf >, <j), <j)', <j>", 
B{(f>,(j),(j)',(j>',</)",t) = 0. Discretizing the operator £ on 
our grid it becomes a non-linear equation for cf>(yi,t n ). 
Our second method is to solve the equation 


C. Time-steps and initial conditions 

The time-step Ax is initially fixed, but can be adjusted 
adaptively depending on how close the denominator B 
in the equation for ^ = A j s t 0 ze ro: we cut the log- 
time-step proportional to the value of the denominator. 
This only applies for the void simulations of the Cubic 
Galileon. In all other cases we used a fixed Ax. The 
initial values we choose depend on the application, but 
in general we just put <f> = 0 and Q = Q at the initial 
time-step. 

Most of our numerical solutions have been derived us¬ 
ing both of the integration methods to cross-check the 
results. In general the NGS solver was found to be the 
most stable one; allowing larger time-steps than the leap¬ 
frog without breaking down, but then again it is much 
much slower. 


D. Treatment of the matter sector 


£(<KVi,t n )) = 0 (52) 


at each time-step using Newton-Gauss-Seidel relaxation 


r ew (yi,tn) = cf> oia (yi,t n )- 


£(Vi, ^n) 


9£{yi,t n )/d<j>{yi,t n ) 


(53) 


with red-black colouring of the grid nodes. Note that 
this method is much slower than our main method (the 
leap-frog) as we must solve an algebraic equation at every 
step as opposed to simply updating a value. 

Time-derivatives, such as <j> and cf>, are calculated using 
the backward-scheme 




l) 

At 

4>(t n ) - 2(j>(t n -i) + </>(t n - 1) 
(At) 2 


(54) 

(55) 


The evolution (gravitational collapse) of matter will 
in general be affected by both the standard Newtonian 
force and the fifth-force from the scalar fielcJ3 Both of 
these terms should be calculated simultaneously. How¬ 
ever, since the purpose of this paper is to study the evo¬ 
lution of the scalar field we have simply chosen to im¬ 
pose a density profile and then study the evolution of 
the scalar field given this profile. This simplifies the nu¬ 
merical solution as we don’t have to evolve the matter 
sector. This is a self-consistent procedure as one can eas¬ 
ily imagine (in principle) that the profile is being set up 
and held together by non-gravitational forces. We choose 
the mexican-hat profile 


5 m {a,r) = S 0 



e ^f(a) 


(58) 


for second order accuracy. Grid derivatives are calculated 
as 


4>"{yi) 


j%»+i) - 4>{yi- 1 ) 

2A y 

Hyi+i) ~ ^ivi) + <t>(yj- 1 ) 

(Ay) 2 


(56) 

(57) 


where /(a) is some function regulating the amplitude of 
the matter perturbation as a function of time. This pro¬ 
file has the advantage that / 0 °° 5 m (t, r)47rr 2 dr = 0 so our 
simulations will not have any excess matter compared to 
the cosmic mean (which could possibly bias the results). 
Unless stated otherwise the main choice for / is /(a) = a 3 
for time-evolving profiles and /(a) = 1 for static profiles. 


7 We have also tried using fifth-order splines, but this is much 
slower and was not found to significantly improve the accuracy 
of the solution. 


Matter moves under the force given by the potential d', but in 
the non-relativistic weak-field limit we can separate the potential 
't into the standard Newtonian part 'h y and a fifth-force oc V 0. 
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FIG. 3. Evolution of the Noether charge, Eq. ( |45| ), for a 
’cluster’ simulation. The spikes happen when outgoing waves 
reflect of the boundary in our simulation box. However, after 
the reflection the charge returns back to the value where it 
started beforehand. 


E. Boundary conditions 


The Lagrangian for our models does not have any 
mass-term which means that, in the linear limit, the 
scalar field will have waves that are only weakly damped 
(by friction terms that depend on the profile we have 
at any moment) as the scalar waves propagate. This 
can cause a problem when waves are created in the box 
and start propagating out towards the boundary of the 
simulation-box. If we employ standard boundary condi¬ 
tions ^ = 0 at r = R the outgoing waves will reflect at 
the boundary and start traveling inwards. 


The reflection of these can be a potential problem and 
there are a few ways to mitigate it. The simplest brute- 
force solution is to make the box so large that the waves 
won’t have time to reflect back in on the time scale of 
the problem we are solving. A more sophisticated solu¬ 
tion is to try to construct non-reflecting boundary condi¬ 
tions or even introduce artificial damping terms into the 
equations of motion. We have opted for the first, and 
simplest, approach and have investigated its validity by 
using the conservation of the Noether current Eq. (45) 
during the evolution. In Fig. © we show the conserved 
charge for a test-case where a mexican-hat profile (see 
the previous section) with <5o = 100. As the outgoing 
waves reflect from the boundary the Noether charge is 
seen to change, but returns to its previous (small) value 


after the reflection has ended and the wave is traveling 
back in again. We have also explicitly checked, by run¬ 
ning simulations using a very large box so that outgoing 
waves never get back in again. We have found that the 
main results of this paper do not depend critically on this 
choice of boundary condition. 



FIG. 4. The solution for EQ — ^ in a simulation where 

8m = 0. This effectively means that we are solving the back¬ 
ground equation of motion in each grid-cell independently. 
We plot the average value (EQ) (symbols) together with the 
analytical result (solid line). 


V. CODE TESTS 

Here we present some of the tests we have performed 
to ensure that our code is working correctly. 


A. No perturbations 

The first test we perform is to put S m = 0. This effec¬ 
tively solves the background equation of motion in every 
grid cell and Q = Q is the expected result. The result 
displayed in Fig. Q shows perfect agreement with ex¬ 
pectations. 


B. A small perturbation 

Building on the previous test we introduce a small 
perturbation in the field and study the propagation of 
the waves it creates. Linearizing the equation of motion 
around the cosmological background we find 
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8<t> + 12 E<2j - a 2 \7 2 y S0 + 8 EQ - 2c 3 Q 4 + 4 H^j 

+50 ( 12 E 2 ^- + + 36 QE 2 + 6 Q^- - c 2 Q 3 - 18£c 3 Q 4 ) = 3 Q 2 a~ 3 n m 6 m 

\ ax C 3 ax J 


(59) 


The corresponding equation(s) for DGP can be found in 
Appendix[B] From this we can read of the speed of sound, 
in physical coordinates, as 


c 2 = % + 8EQ-2c 3 Q +4 E<§ 


% + 12 EQ 


(60) 


for the Cubic Galileon. We test the code by taking 
6 m = 0 , placing the system in the cosmological back¬ 
ground attractor Q = Q and then add a small, one- 
wavelength, perturbation 50 oc sin(27r y ^* )0(y — y *), 
where 0 is the Heaviside function. 

Solving the wave equation analytically is not easy as 
both the speed of sound and the friction term depends 
on time. However, we can derive an analytical approxi¬ 
mation, 


/ 


E(x) 


sin [-r(y ~ y* ± ct(x)) 

, V dy 

x0(y ~ V* ± ct{x))6{y * + dy + cf(x) - y) (61) 


where ct(x) = J e(x) anc ^ D{x) the t erm i n front of 

S(j) divided by the term in front of S<j>. 

Since the approximation above is not exact we don’t 
expect perfect agreement. We therefore also choose to 
implement and solve the linear field equation Eq. (59) 
and use this to compare our result with (using both of 
our integration methods). In Fig. ^ we show the result 
of this test. The full equations agrees perfectly with the 
linearized solution and also fairly well with our analytical 
approximation. 


grid the further out the wave will start (the minimum 
r/R we can represent in the code is (r/R ) m ; n = j^) thus 
giving it a head-start when comparing it with the more 
refined simulations. 

Next we have tested how the results depend on the 
size of the time-steps used in the computation. First 
of all we should mention that we are required to take 
very small time-steps in order for the solver not to break 
down. The leapfrog solver is more sensitive to the size 
of the time-steps than the NGS solver. The tests shown 
here are for the NGS integrator which allows us to go all 
the way down to N = 100 time steps between a = 0.1 and 
a = 1.0. The results are shown in Fig. 0 - With very few 
time-steps we see a difference in the solution compared to 
the simulations with the most time-steps. This difference 
seems to go away with time and the reason for this is that 
the quasi-static solution is an attractor for the system, 
however we are unable to resolve the oscillations of the 
field if we keep the number of time-steps small (small 
here means less than 1000). As we increase the number of 
time-steps the solution we find stops changing (including 
the oscillations we see in the field). For this particular 
test we found that roughly 10 5 time-steps is needed to 
accurately track the evolution of the field at all times 
from a = 0.1 till a = 1 . 0 . 

Finally we tested the difference between the two dif¬ 
ferent integration methods keeping the time-steps equal. 
In Fig. ([ 8 ]) we show the difference in the solution for <f>{r) 
between the leapfrog and the NGS solver for simulations 
of a growing void. The agreement is excellent where it 
should be. It is only in regions where </> « 0 that we see 
a small difference. 


C. Convergence of the solution 

We have performed several tests to check how sensitive 
our results are to the grid size, the size of the time-steps 
and the integration method used in the calculation. The 
results shown below are simulations of a growing mexican 
hat void profile (/(a) = a 3 ) with final density contrast 
<5 0 = -0.5. 

In Fig. (Joj) we show results from calculations using N = 
512,1024, 2048 and 4096 grid nodes respectively. The 
results agree very well across the different resolutions. 
The only real difference we found is in the first few time- 
steps where scalar waves are emitted (since we start off 
with (f> = 0 ) going out from the object at early times. 
The waves are created at r/R = 0 and the coarser the 


VI. SIMULATION RESULTS 

In this section we will present the results of our simu¬ 
lations and try to answer the questions we posed in the 
introduction. 


A. Does the Vainshtain solution always forms? 

To try to answer the question in the title we impose 
an analytic density profile S m (x, y) and release the scalar 
field at <fi = 0 with Q = Q and follow the subsequent 
evolution. 

In all cases we looked at, for both over- and under¬ 
densities, the profile quickly (with the speed of sound) 
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FIG. 5. The evolution of a one-wavelength sine wave perturbation 5(j>, starting from r/R = 0.5. We show both the ingoing 
and outgoing wave, at three snapshots (in the order blue, green and then red) after the initial release together with our 
analytical approximation Eq. (611 for the self-accelerating DGP model with r c Ho = 1.35. The difference between the full 
non-linear solution and the solution found by implementing the linearized evolution equation is indistinguishable with both of 
our methods of integration. 


evolved towards the quasi-static solution and started os¬ 
cillating around it until it finally settled down. In Fig. Q 
we show the evolution as a function of time for the case of 
a growing void with density contrast do = —0.5ata = 1.0 
and in Fig. (101 we show the evolution for a growing clus¬ 
ter with density contrast do = 100 at a = 1.0. 

We have also tried experiments with other initial con¬ 
ditions with the same result. The only case where we 
saw a real difference in the full solution compared to the 
quasi-static analytical approximation was when we con¬ 
sidered the case of an object several giga-parsec large 
growing very rapidly in time. In that, unphysical, case 
the speed of sound of the field is not large enough for the 
field to be able to evolve quickly enough to catch up with 
the evolving quasi-static solution. 


B. Is the Vainshtein solution stable? 

To test the stability of the Vainshtein solution we run 
simulations where we have a fixed density profile where 
the scalar field has relaxed at the quasi-static solution 
and then we send in waves in the scalar field traveling 
towards the object. To form these waves we take advan¬ 
tage of reflecting boundary conditions to get the waves 
created when the scalar field evolved from <fi = 0 and 
into the Vainshtein solution. By using a large enough 
box we can get these waves to reflect from the boundary 
and return to the center at any give time we want. 

When scalar waves hit the object they are quickly re¬ 
flected back out again. The quasi-static solution remains 
stable under this process, changing very little, and only 


when we hit the profile with waves of very large ampli¬ 
tude are we able to significantly change the profile for a 
little while. However after the waves have been reflected 
the field profile relaxes to its old position, see Fig. (13). 
The Vainshtein solution seems very stable to interference 
from incoming scalar waves. 

It would be interesting to redo this experiment for the 
case where the matter making up the density profile actu¬ 
ally experiences the fifth-force present and study if scalar 
waves can have any impact on the density distribution. 
This is however beyond the reach of this paper. 


C. How good is the quasi-static approximation? 


To start with the punchline: the quasi-static approxi¬ 
mation is excellent in almost all cases we have checked. 
If we start the simulation with </> being away from the 
quasi-static solution then it will quickly start to evolve 
towards it and once it has been reached, it undergoes 
damped oscillations until it quickly settles down to it 
(as shown in the previous section). If the quasi-static 
solution is evolving (like it will be if our density pro¬ 
file is growing/decaying with time) then we find that the 
quasi-static solution acts as an attractor which the full 
numerical solution follows closely. 

In Fig. © and Fig. ( [To| (see the lower left panel) we 
show the evolution of w together with the quasi-static 
solution for the ubic Galileon for a void and a cluster 
respectively. The same figure for the DGP model can be 
seen in Fig. (121 and Fig. (111. 


The deviations from quasi-static evolution we find for 
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FIG. 6. Convergence test with respect to the number of grid nodes we use. Here we show the profiles <p(y) (top), SQ = H 0 1 
(middle) and ui = R ^ at a = 1.0 (left) and a = 0.1 (left) for four simulations using N = 512,1024, 2048 and 4096 grid nodes 
respectively. The results agree very well at both times and all radii, with only some small differences visible for the N = 512 
case compared to the more refined simulations. These results are for the Cubic Galileon, but we get the same behavior for the 
DGP model. 


the DGP model are consistent with the simulations of 
[22] where the effect of the approximation was estimated 
from the time-evolving quasi-static solution (but with¬ 
out actually having any time evolution in the simulation 
itself). 

It would be interesting to study what happens when 
the density profile changes dramatically in a short period 
of time, but this is beyond what we aim to study in this 


paper. 


D. Validity of other approximations 

We have run simulations where we evolve the field 
equation both with and without taking the metric deriva¬ 
tive terms 4>, \1 >,d$/dr and /dr into account in the 
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0.0001 0.001 0.01 0.1 1 
r/R (R = 500 Mpc/h) 


FIG. 7. Convergence test with respect to the number of time-steps we use. Here we show the held profile (f>(r ) at a = 0.15,0.2 
and a = 0.4 for five different simulations using N t i me = 100, 500,1000, 5000 and 10000 time-steps between a = 0.1 and a = 1.0. 
At early times, when the held is still evolving towards the quasi-static solution, the error can be quite large if we use too few 
time-steps. However, at later times when the held have settled close to the (evolving) quasi-static solution this difference has 
largely been washed away. These results are for the Cubic Galileon, but we get the same behavior for the DGP model. 


field equation. Note that neglecting terms like <f>, ^ are 
strictly speaking part of the usual quasi-static approxi¬ 
mation, but we have chosen to study it separately here. 

In the equation of motion we do enforce this approxi¬ 
mation by simply setting all the A and B terms in Eq. (36 
and Eq. (B2) to zero. The difference in the fifth-force 


E. Breakdown of the numerical solution in deep 
voids 


w = R(ir between the two simulations, with and without 
this approximation, is shown in Fig. @. 

The difference in the force inside the overdensity is 
seen to be at the 10 -5 level which is of the same order 
of magnitude as the Newtonian potential IF of the ob¬ 
ject in agreement with the rough arguments we gave in 
Sec. (III C). Outside the overdensity, on the other hand, 
the difference can be as large as 10 -2 . However we should 
note that in the region (r > 10 Mpc/h from the center 
of the overdensity) where we find the largest difference, 
the force itself is very close to zero, see Fig. (10 1 , so the 


physical effect this difference would have on matter lo¬ 
cated in this region is practically as small as inside the 
object. For the underdensity we simulated we get a very 
similar result. 


Lastly we turn to the case discussed in Sec. ( |HDt . 
For the Cubic Galileon, at late times and in deep voids 
the quasi-static approximation ceases to exist, i.e. we get 
complex solutions in Eg. (25 ). When we reach the regime 


where the quasi-static approximation breaks down, so do 
our simulations. In particular, we find that the denomi¬ 
nator B in the field equation = 5 crosses zero, lead¬ 
ing to nonsensical solutions. If B —)■ 0 then Q —> 00 


and the speed of sound also diverges, see Eg. (36 1 , un¬ 
less A —> 0 at the same time. In practice we stopped 
the simulation if B crossed zero, and if we continue the 
simulation past this point then the profile quickly de¬ 
velops kinks and breaks down (in this regime the sound 
speed squared is negative). This effect happens in our 
simulations no matter how small we take the time-steps 
to be. We have gone down to time-steps of the order 
Sa ~ 10 -10 corresponding to 1 week (!) of cosmic time 
per time-step and still the solution crosses zero without 
any sign of turning around. This problem is encountered 
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FIG. 8. The solution for the Cubic Galileon (f)(r ) normalized 
such that <j>( 0) = 1 together with the difference between the 
solutions found by using the two different integration methods 
(leapfrog and NGS). For this test we used a IV = 4096 grid 
with IVtime = 5 ■ 10 4 time-steps from a = 0.1 till a = 1.0 
and where the density profile was a growing void with 8 o = 
— 0.5. We see that the difference stays below the ~ 10 -4 level 
whenever <j> itself is larger than ~ 1CF 4 times its maximum 
value. 


for both methods of integration and with all gridsizes, 
from N = 100 to N = 10000, we have tried. 


There seems to be no terms in the equation of mo¬ 
tion that are able to stop B from crossing zero and even 
though one should be careful in extracting definite con¬ 
clusions from a failed numerical experiment we think it 
is, at least, very unlikely that the quasi-static approxima¬ 
tion is able to alleviate this problem and that it is indeed 
a true instabillity of the Cubic Galileon model. 


From a more theoretical point of view, a similar type 
of instabillity was discussed in m for the G-bounce sce¬ 
nario. 


For the DGP model we do not have this problem and 
both the numerical solution and the quasi-static solution 
exists for all voids at all times. 


VII. SUMMARY 

We have studied the Vainshtein mechanism without 
restricting ourselves to the quasi-static approximation. 
By solving for the full time-evolution of the DGP and 
the Cubic Galileon scalar field in a spherical symmetric 
spacetime we are able to demonstrate that the quasi¬ 
static approximation is indeed a very good approxima¬ 
tion. This result has previously been explicitly shown 
for several other modified gravity models and with this 
study there are now explicit checks of this approxima¬ 
tion for all the major screening mechanisms present in 
the literature. 

We found that the quasi-static Vainshtein solution is a 
stable attractor for the evolution of the scalar field. It is 
an attractor in the sense that even if we release the scalar 
field far away from this solution it will quickly evolve into 
and relax to this solution. It is stable in the sense that it 
is hard to disrupt it: incoming scalar waves are absorbed 
by the profile and quickly emitted back again without 
significantly changing the profile. 

We have also demonstrated, as expected, that other 
approximations closely linked to the quasi-static approx¬ 
imation are also valid. For example we found that ne¬ 
glecting terms proportional to time and single spatial 
derivatives of the metric potentials (<h, d<&/dr) in the field 
equation for </> have very little impact on the solution. 

The final question we set out to answer was if relaxing 
the quasi-static approximation could aleviate the prob¬ 
lem found in |2T. regarding the breakdown of the solution 
in deep voids close to the present time. By solving for the 
full evolution of the scalar field we encountered exactly 
the same problem. Our numerical solution breaks down 
at roughly the same point in time independent of the 
time-steps, gridsize and numerical integrator we choose. 
It therefore seems very likely that we are dealing with a 
true instabillity of the model. 
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Appendix A: The A terms 

Here we show the expressions for the Ai-q factors in 
the equation of motion Eq. (36 ): 
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FIG. 9. The evolution of the Cubic Galileon field (ft an< i its derivatives from a = 0.1 to a = 1.0 for a void with under-density 
Jo = —0.5 today. The dashed lines in the to plot shows the analytical quasi-static solution. 
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Appendix B: Equations for the DGP model 


Here we present, for completeness, the DGP equations 
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FIG. 10. The evolution of the Cubic Galileon scalar field <f> an< l its derivatives from a = 0.1 to a = 1.0 for a cluster with over 
density So = 100 today. The dashed lines in the uj plot shows the quasi-static solution. 




FIG. 11. The gradient of <f>, uj = at three different times for a cluster with So = 100 at the present time ( S m (r, a) oc a 3 ) for 
the normal branch (left) and the self-accelerating branch (right) of the DGP model. In both cases we have taken rcHo = 1.35 
which in the self-accelerating branch corresponds to = 0.26. The solid lines shows the analytic quasi-static solution. 
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FIG. 12. The gradient of <j>, u> = at three different times for a void with So = —1.0 at the present time ( S m (r,a ) oc a 3 ) for 
the normal branch (left) and the self-accelerating branch (right) of the DGP model. In both cases we have taken rc-ffo = 1.35 
which in the self-accelerating branch corresponds to fl m = 0.26. The solid lines shows the analytic quasi-static solution. 
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In the DGP model the background value for the scalar 
field is simply <5 = 0. The Hubble equation for the self- 
accelerating branch is given by 


E(a) = 


1 


2 r c H 0 


+ A / H m d 3 + 


1 


(2 r c H 0 ) 


(B6) 


The linearized equation of motion reads 


t dSQ 
dx 


E(a) + 3 E(a)SQ - = 0 (B8) 


a 2 /3(a) 


where <50=^^- 


where Q m = 1 — r ^ or equivalently t c Hq = . 

For the normal branch we assumed dark-energy in such 
a form to give us a ACDM background 


E(a) = v fi m a 3 + Da 
where Ha = 1 — H m . 


Appendix C: Useful equations 


In this appendix we present some very useful equations 
(B7) needed to derive the field equation for a scalar field theory 
for spherical symmetry and in a perturbed, to first order 
in perturbation theory for the metric potentials, FRLW 
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FIG. 13. We show the response of the ui profile to a scalar wave (coming in from the right). The wave reflects from the profile 
which then quickly returns to the quasi-static solution. 


metric in the Newtonian gauge: 

ds 2 = — (1 + 2^/)dt 2 + a 2 (l — 2<I>) x 

(dr 2 + r 2 dd 2 + r 2 sin 2 9d<j) 2 ) (Cl) 

Below we use the notation {1,2,3,4} = {t, r, 9 , </>}, a dot 
is a time-derivative and a prime is a radial derivative 
' = d/dr. To simplify the notation we use V 2 for the 
radial Laplace operator: V 2 </> = </>" + 2 cj>'. 

We have the following relations for the Christoffel sym- 


The T 1 terms reads 


r 1 

1 ii 

= 4- 

(C2) 

r 1 

1 12 

= T' 

(C3) 

r 1 

1 22 

= Ha 2 - 2 Ha 2 (<f> + \E') - a 2 $ 

(C4) 

r 1 

1 33 

= r 2 a 2 H - 2 r 2 a 2 H(<5> + T) - r 2 a 2 $ 

(C5) 

r 1 

1 44 

= r 2 a 2 H sin 2 9 - r 2 a 2 sin 2 9(AH + T-) 

(C6) 
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FIG. 14. The difference in the force u> = for a simulation where we included metric derivatives (4>, d<f?/<7r etc.) compared 
to one where we neglected these. The simulations are of a cluster with overdensity <5o = 100 (left) and a void of underdensity 
Jo = —0.5 (right) at the present time. 


The T 2 terms reads 


2 4-' 

111 " a 2 

(C7) 

r ?2 = 77 - $ 

(C8) 

r 22 = 

(C9) 

r 33 = -r + r 2 4-' 

(CIO) 

T 2 4 = -r sin 2 8 - 2r sin 2 0(4- — 4-) + r 2 sin 2 04-' 

(Cll) 

The T 3 terms reads 


r?3 = H - 4> 

(C12) 

ri 3 = 

(C13) 


(C14) 

The T 4 terms reads 


F[ 4 = 77 - 4- 

(CIS) 

r 24 = ' - *' 

(C16) 

Tg 4 = COt 0 

(C17) 


The rest of the non-zero terms follows from using the 
symmetry T“„ = on the terms presented above. 

The determinant of the metric g is 

det g = —r 4 a 6 sin 2 9(1 — 4$) (C18) 


Some useful components of the Ricci tensor are 


Ru = 3- - - 4774- - 5774- - 2$ - 4- (C19) 

a a 2, 

R 22 = -2a 2 77 2 - a 2 -( 1 - 2$ - 24-) + 4a 2 77 2 (4» + 40 
a 

- V 2 4> + 5a 2 77$ + 2a 2 774- + a 2 $ (C20) 


R 12 = -2 774-' 


4> - 4- 


— 4 ,/ — 4- / 


(C21) 


The Ricci scalar becomes 


R = - 6 


77 2 


(1 - 24') - -^V 2 4- 


+ -($'- 4'') + 16774- + 14774- + 44- + 24- 


ra 


(C22) 


Some useful derivatives of the scalar field is 

= ^(l + 2$)-^ 2 (l-24-) (C23) 

□(/) = (1 + 9 2$) V 2 </> - 3770(1 - 24-) - 4> (C24) 
a z 


r 
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